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Abstract 

A novel Stochastic Event-Driven Molecular Dynamics (SEDMD) algorithm is developed for the 
simulation of polymer chains suspended in a solvent. The polymers are represented as chains 
of hard spheres tethered by square wells and interact with the solvent particles with hard core 
potentials. The algorithm uses Event-Driven Molecular Dynamics (EDMD) for the simulation of the 
polymer chain and the interactions between the chain beads and the surrounding solvent particles. 
The interactions between the solvent particles themselves are not treated deterministically as in 
event-driven algorithms, rather, the momentum and energy exchange in the solvent is determined 
stochastically using the Direct Simulation Monte Carlo (DSMC) method. The coupling between 
the solvent and the solute is consistently represented at the particle level, however, unlike full 
MD simulations of both the solvent and the solute, the spatial structure of the solvent is ignored. 
The algorithm is described in detail and applied to the study of the dynamics of a polymer chain 
tethered to a hard wall subjected to uniform shear. The algorithm closely reproduces full MD 
simulations with two orders of magnitude greater efficiency. Results do not confirm the existence 
of periodic (cycling) motion of the polymer chain. 



I. INTRODUCTION 



Driven by nanoscience interests it has become necessary to develop tools for hydrody- 
namic calculations at the atomistic scale [H El [3l [U [5] . Of particular interest is the modeling 
of polymers in a flowing "good" solvent for both biological (e.g., cell membranes) and en- 
gineering (e.g., micro-channel DNA arrays) applications [H [6j. The most widely studied 
polymer models are simple linear bead-spring; freely-jointed rods; or worm-like chains. Such 
models have been parameterized for important biological and synthetic polymers. Much 
theoretical, computational, and experimental knowledge about the behavior of these models 
has been accumulated for various representations of the solvent. However, the multi-scale 
nature of the problem for both time and length is still a challenge for simulation of reason- 
ably large systems over reasonably long times. Furthermore, the omission in these models 
of the explicit coupling between the solvent and the polymer chain (s) requires the introduc- 
tion of adjustable parameters (e.g., friction coeflicients) to be determined empirically. The 
algorithm presented here overcomes this for a linear polymer chain tethered to a hard wall 
and subjected to a simple linear shear flow [3 [U [9l [IHl [H]. Of particular interest is the 
long-time dynamics of the polymer chain |7l [9l [IHl [121 113] and any eflFects of the polymer 
motion on the flow fleld. 

Brownian dynamics is one of the standard methods for coupling the polymer chains to 
the solvent [HI [15] . The solvent is only implicitly represented by a coupling between the 
polymer beads and the solvent in the form of stochastic (white-noise) forcing and linear 
frictional damping. The flow in the solvent is not explicitly simulated, but approximated 
as a small perturbation based on the Oseen tensor. This approximation is only accurate at 
large separations of the beads and at sufliciently small Reynold's numbers. Even algorithms 
that do model the solvent explicitly via Lattice Boltzmann (LB) [16j, incompressible (low 
Reynolds number) CFD solvers [TTl [181 [iH] , or multiparticle collision dynamics [2Ql [2U [22l 
[23] , typically involve phenomenological coupling between the polymer chain and the flowing 
fluid in the form of a linear friction term based on an effective viscosity (a notable exception 
being the algorithm described in Ref. [23]). Furthermore, solvent fluctuations in the force on 
the polymer beads are often approximated without fully accounting for spatial and temporal 
correlations. Finally, the reverse coupling of the effect of the bead motion on the fluid flow 
is either neglected or approximated with delta function forcing terms in the continuum fluid 
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solver [24] . More fundamentally, continuum descriptions of flow at micro and nanoscales are 
known to have important deficiencies [K and therefore it is important to develop an all- 
particle algorithm that is able to reach the long times necessary for quantitative evaluation 
of approximate, but faster, algorithms. 

The most detailed (and expensive) modeling of polymers in flow is explicit molecular dy- 
namics (MD) simulation of both the polymer (solute) and the surrounding solvent pTl [25] . 
Multi-scale algorithms have been developed to couple the MD simulation to Navier-Stokes- 
based computational fluid dynamics (CFD) calculations of the flow fleld [8]. However, the 
calculation time still remains limited by the slow molecular dynamics component. Thus 
the computational eflFort is wasted on simulating the structure and dynamics of the solvent 
particles, even though it is the polymer structure and dynamics (and their coupling to the 
fluid flow) that is of interest. Our algorithm replaces the deterministic treatment of the 
solvent-solvent interactions with a stochastic momentum exchange operation, thus signifl- 
cantly lowering the computational cost of the algorithm, while preserving microscopic details 
in the solvent-solute coupling. 

Fluctuations drive the polymer motion and must be accurately represented in any model. 
Considerable effort has been invested in recent years in including fluctuations directly into 
the Navier-Stokes (NS) equations and the associated CFD solvers [5l[I7l[2S|. Such fluctuat- 
ing hydrodynamics has been coupled to molecular dynamics simulations of polymer chains 
[l9], but with empirical coupling between the beads and the fluid as discussed above. To 
avoid the empirical coupling, the solvent region could be enlarged by embedding the atom- 
istic simulations of the region around the polymer chain (such as pure MD or our combined 
MD/DSMC algorithm) in a fluctuating hydrodynamics region. The bidirectional coupling 
between the continuum and particle regions has to be constructed with great care so that 
both fluxes and fluctuations are preserved [27j. A well-known problem with such multiscale 
approaches is that the flnest scale (atomistic simulation) can take up the majority of com- 
putational time and thus slow down the whole simulation. By using DSMC the cost of the 
particle region can be made comparable to that of the continuum component. 

The Stochastic Event-Driven Molecular Dynamics (SEDMD) algorithm presented here 
combines Event-Driven Molecular Dynamics (EDMD) for the polymer particles with Direct 
Simulation Monte Carlo (DSMC) for the solvent particles. The polymers are represented as 
chains of hard spheres tethered by square wells. The solvent particles are realistically smaller 
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than the beads and are considered as hard spheres that interact with the polymer beads with 
the usual hard-core repulsion. The algorithm processes true (deterministic, exact) binary 
collisions between the solvent particles and the beads, without any approximate coupling or 
stochastic forcing. However, for the purposes of the EDMD algorithm, the solvent particles 
themselves do not directly interact with each other, that is, they can freely pass through each 
other as for an ideal gas. Deterministic collisions between the solvent particles are replaced 
with stochastic DSMC collisions. Both asynchronous (event-driven) and synchronous (time- 
driven) algorithmic ways of processing these stochastic collisions will be discussed in the 
next section. Note that our algorithm is similar to a recent algorithm developed for soft 
interaction potentials combining time-driven MD with multiparticle collision dynamics [23j . 

The fundamental ideas behind our algorithm are described next, and further details are 
given in Section |III[ Section M gives results from the application of the algorithm to the 



tethered polymer problem, and some concluding remarks are given in Section VI 



II. HYBRID COMPONENTS 



In this section we briefly describe the two components of the SEDMD algorithm: The 
stochastic handling of the solvent and the deterministic handling of the solute particles. 
These two components are integrated (i.e., tightly coupled) into a single event-driven algo- 
rithm in Section lllli 



A. Solvent DSMC Model 



The validity of the (incompressible) Navier-Stokes continuum equations for modeling mi- 
croscopic flows has been well established down to length scales of 10 — lOOnm [3j. However, 
there are several issues present in microscopic flows that are diflicult to account for in mod- 
els relying on a purely PDE approximation. Firstly, it is not a priori obvious how to treat 
boundaries and interfaces well so as account for the non-trivial (possibly non-linear) cou- 
pling between the flow and the microgeometry. Furthermore, fluctuations are not typically 
considered in Navier-Stokes solvers, and they can be very important at instabilities [28] or 
in driving polymer dynamics. Finally, since the grid cell sizes needed to resolve complex 
microscopic flows are small, a large computational eflFort (comparable to DSMC) is needed 
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even for continuum solvers. An alternative is to use particle-based methods, which are ex- 
plicit and unconditionally stable and rather simple to implement. The solvent particles are 
directly coupled to the microgeometry, for example, they directly interact with the beads 
of a polymer. Fluctuations occur naturally with the correct spatio-temporal correlations. 
However, as in continuum descriptions, the structure of the fluid is lost and under certain 
conditions the high compressibility of the DSMC (ideal gas) fluid can cause difliculties. 

Several particle methods have been described in the literature, such as MD p5], dissipative 
particle dynamics (DPD) [29], and multi-particle collision dynamics (MPCD) [2l[23j. Our 
method is similar to MPCD (also called stochastic rotation dynamics or the Malevanets- 
Kapral method), and in fact, both are closely related to the Direct Simulation Monte Carlo 
(DSMC) algorithm of Bird [30j. The key idea behind DSMC is to replace deterministic 
interactions between the particles with stochastic momentum exchange (collisions) between 
nearby particles. Speciflcally, particles are propagated by a flxed time step At, as in MD, 
moving ballistically along straight lines during a time-step (advection step). At the end 
of each time step, the particles are sorted into cells, each containing on the order of ten 
particles, and then a certain number of random pairs of particles that are in the same cell 
are chosen to undergo stochastic collisions (collision step). These collisions do not take into 
account the positions of the particles other than the fact they are in the same cell (i.e., they 
are nearby). The collisions conserve momentum and energy (but not angular momentum) 
exactly. Formally, DSMC can be seen as a method for solving the Boltzmann transport 
equation for a low-density gas, however, it is not limited to gas flows [HU [32l [33]. Our 
purpose for using DSMC is as a replacement for expensive MD, preserving the essential 
hydrodynamic "solvent" properties: local momentum conservation, and linear momentum 
exchange on length scales comparable to the particle size, and a similar fluctuation spectrum. 

In the multiparticle collision variant of this algorithm originally proposed by Kapral, the 
traditional DSMC collection of binary collisions is replaced by a multi-particle collision in 
which the velocities of all particles in the cell are rotated by a random amount around the 
average velocity [2, 23j. This change improves efliciency but at the cost of some artiflcial 
effects such as loss of Galilean invariance. These problems can be corrected and the method 
has been successfully used in modeling polymers in flow by including the beads, considered 
as (massive) point particles, in the stochastic momentum exchange step [20l |22l |34]. We 
will employ traditional DSMC in our algorithm in order to mimic the actual (deterministic) 
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momentum exchange between solvent molecules (as it would be in an MD simulation) and 
in order to avoid any possible artifacts. 

A fundamental deficiency of DSMC as a (micro or nano) hydrodynamic solver is the 
large (ideal gas) compressibility of the fiuid. For subsonic fiows this compressibility does not 
qualitatively affect the results as the DSMC fiuid will behave similarly to an incompressible 
liquid, however, the (Poisson) density fiuctuations in DSMC are significantly larger than 
those in realistic liquids. Furthermore, the speed of sound is small (comparable to the average 
speed of the particles) and thus subsonic (Mach number less than one) fiows are limited to 
relatively small Reynolds numbers ^. The Consistent Boltzmann Algorithm (CBA) [32l [33] . 
as well as algorithms based on the Enskog equation [35l[36], have demonstrated that DSMC 
fiuids can have dense-fiuid compressibility. A similar algorithm was recently constructed for 
MPCD [37j . We are currently evaluating several DSMC variants in terms of their efiiciency 
and thermodynamic consistency under high densities ^ and will report our findings in future 
work. 

B. Polymer MD Model 

Polymer chains in a solvent are modeled using continuous pair potentials and time-driven 
MD (TDMD), in which particles are synchronously propagated using a time step At, inte- 
grating the equations of motion along the way. For good solvents, the polymer beads are 
represented as spherical particles that interact with other beads and solvent particles with 
(mostly) repulsive pair potentials, such as the positive part of the Lennard- Jones potential. 
Additionally, beads are connected via (usually finitely-extensible FENE or worm-like) springs 
in order to mimic chain connectivity and elasticity [25j . Additionally, stochastic forces may 
be present to represent the solvent. The time steps required for integration of the equations 
of motion in the presence of the strongly repulsive forces is small and TDMD cannot reach 
long time scales even after parallelization. An alternative is to use hard spheres instead 

^ For a low-density gas the Reynolds number is Re = M/K^ where M = Vfiow/c is the Mach number, and 
the Knudsen number K = X/L is the ratio between the mean free path A and the typical obstacle length 
L. This shows that subsonic flows can only achieve high Re flows for small Knudsen numbers, i.e., large 
numbers of DSMC particles. 

^ Note that the density fluctuations in the CBA fluid are identical to those in an ideal gas and thus 
thermodynamically inconsistent with the compressibility. 
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of soft particles, allowing replacement of the FENE springs with square- well tethers, thus 
avoiding the costly force evaluations in traditional MD. Hard sphere MD is most efficiently 
performed using event-driven molecular dynamics (EDMD) [38l |39l HOI [41]. If the detailed 
structure and energetics of the liquid is not crucial, such EDMD algorithms can be just as 
effective as TDMD ones but considerably faster. The essential difference between EDMD 
and TDMD is that EDMD is asynchronous and there is no time step, instead, collisions 
between hard particles are explicitly predicted and processed at their exact (to numerical 
precision) time of occurrence. Since particles move along simple trajectories (straight lines) 
between collisions, the algorithm does not waste any time simulating motion in between 
events (collisions). 

Hard-sphere models of polymer chains have been used in EDMD simulations for some 
time [lOl nil [43j . These models typically involve, in addition to the usual hard-core exclusion, 
additional square well interactions to model chain connectivity. The original work by Alder 
et al on EDMD developed the collisional rules needed to handle arbitrary square wells [38j . 
Infinitely high wells can model tethers between beads, and the tethers can be allowed to be 
broken by making the square wells of finite height, modeling soft short-range attractions. 
Recent studies have used square well attraction to model the effect of solvent quality [41j. 
Even more complex square well models have been developed for polymers with chemical 
structure and it has been demonstrated that such models, despite their apparent simplicity, 
can successfully reproduce the complex packing structures found in polymer aggregation 
[421 [43] . Recent work on coupling a Kramer bead-rod polymer to a NS solver has found 
that the use of hard rods (instead of soft interactions) not only rigorously prevents rod-rod 
crossing but also achieves a larger time step, comparable to the time step of the continuum 
solver [21j. 

This study is focused on the simplest model of a polymer chain, namely, a linear chain 
of Ni) particles tethered by unbreakable bonds. This is similar to the commonly-used freely 
jointed bead-spring FENE model model used in time-driven MD. The length of the tethers 
has been chosen to be where D5 is the diameter of the beads ^. The implementation 

of square- well potentials is based on the use of near-neighbor lists (NNLs) in EDMD, and 

^ Note that the hard-sphere model rigorously prevents chain crossing if the tether length is less than \f2D\) 
since two tethers shorter than this length cannot pass through each other without violating impenetrability. 
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allows for the specification of square- well interactions for arbitrary pairs of near neighbors. 
In particular, one can specify a minimal L^'^^ > D^j and maximal distance (tether length) 
j^max ^ j^mm f^^. arbitrary pairs of near neighbors ^. 

III. DETAILS OF HYBRID ALGORITHM 

In this section the hybrid EDMD/DSMC algorithm, which we name Stochastic EDMD 
(SEDMD), is described in detail. Only a brief review of the basic features of EDMD is 
given and the focus is on the DSMC component of the algorithm and the associated changes 
to the EDMD algorithm described in detail in Ref. [39j. A more general description of 
asynchronous event-driven particle algorithms is given in Ref. [l^- 

Asynchronous event-driven (AED) algorithms process a sequence of events (e.g., colli- 
sions) in order of increasing event time tg. The time of occurrence of events is predicted 
and the event is scheduled to occur by placing it an event queue. The simulation iteratively 
processes the event at the head of the event queue, possibly scheduling new events or inval- 
idating old events. One impending event per particle i, 1 < i < A^, is scheduled to occur at 
time with partner p (e.g., another particle j). The particle position and velocity are 
only updated when an event involving particle i is processed and the time of last update ti is 
recorded (we will refer to this procedure as a particle update). We note that traditional syn- 
chronous time-driven (STD) algorithms with a time step At are a trivial variant of the more 
general AED class. In particular, in an STD algorithm events occur at equispaced times 
and each event is a time step requiring an update of all of the particles. The AED algorithm 
processes a mixture of events involving single particles or pairs of particles with time steps 
that involve the simultaneous (synchronous) update of a large collection of particles. 

Every particle i belongs to a certain specie Si. Particles with species Si and Sj may or may 
not interact with each other (i.e., they may not be subject to the hard-particle non-overlap 
condition). We focus on a system in which a large fraction of the particles belong to a 
special specie Sdsmc representing DSMC particles (e.g., solvent molecules). These DSMC 
particles do not interact with each other (i.e., they freely pass through each other), but 
they do interact with particles of other species. We focus on the case when the non-DSMC 

^ A value LJ^*^ > can be used to emulate chain rigidity (i.e., a finite persistence length) by using second 
nearest-neighbor interactions between chain beads. 
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particles are localized in a fraction of the simulation volume, while the rest of the volume 
is filled with DSMC particles. This will enable us to treat the majority of DSMC particles 
sufficiently far away from non-DSMC particles more efficiently than those that may collide 
with non-DSMC particles. 

Before describing the SEDMD algorithm in detail, we discuss the important issue of 
efficiently searching for nearby pairs of particles. 

A. Near Neighbor Searches 

When predicting the impending event of a given particle i, the time of potential collision 
between the particle and each of its neighbors (nearby particles) is predicted [39l |44] . The 
DSMC algorithm also requires defining neighbor particles, that is, particles that may col- 
lide stochastically during the DSMC collision step. For efficiency, geometric techniques are 
needed to make the number of neighbors of a given particle 0(1) instead of 0{N). 

In SEDMD we use the so-called linked list cell (LLC) method for neighbor searching in 
both the MD and DSMC components. The simulation domain is partitioned into Nceiis cells 
as close to cubical as possible. Each particle i stores the cell q to which its centroid belongs, 
and each cell c stores a list Cc of all the particles it contains, as well as the total number 
of particles Nc in the cell. For a given interaction range, neighbors are found by traversing 
the lists of as many neighboring cells as necessary to ensure that all particles within that 
interaction range are covered. In traditional DSMC, only particles within the same cell are 
considered neighbors and thus candidates for collision. There are also variants of DSMC in 
which particles in nearby cells are included in order to achieve a non-ideal equation of state 
[35l|36j. A more general implementation would use different cell meshes for MD and DSMC 
neighbor searches, however, that would significantly complicate the implementation. 

1. Cell Bitmasks 

In addition to the list of particles each cell c stores a bitmask Aic consisting of 
^bits > + 4 bits (bitfields). These bits may be one (set) or zero (not set) to indicate 
certain properties of the cell, specifically, what species of particles the cell contains, whether 
the cell is event or time driven, and to specify boundary conditions. In order to distinguish 
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the cells that contain non-DSMC particles (i.e., particles of specie other than sdsmc) from 
those that contain only DSMC particles, bit 7 is set if the cell may contain a particle of 
specie 7. The bit is set whenever a particle of specie 7 is added to the cell, and all of the 
masks are reset and then re-built (i.e., refreshed) periodically. When performing a neighbor 
search for a particle i, cells not containing particles of species that interact with specie Si are 
easily found (by OR'ing the cell masks with a specie mask) and are simply skipped. This 
speeds up the processing of DSMC particles since cells containing only DSMC particles will 
be skipped without traversing their lists of particles. 

For the purposes of the combined MD/DSMC algorithm we will also need to distinguish 
those cells that are nearby non-DSMC particles, that is, that contain particles within the 
interaction range of some non-DSMC particle. Such cells will be treated using a fully event- 
driven (ED) scheme, while the remaining cells will be treated using a time-driven or mixed 
approach. We use one of the bits in the bitmasks, bit 7^d, to mark event-driven (ED) 
cells whenever a neighbor search is performed for a non-DSMC particle. Specifically, bit 
"^ED is set for a given cell whenever the cell is traversed during a neighbor search for a non- 
DSMC particle. This scheme correctly masks the cells by only modifying the neighbor search 
routines without changing the rest of the algorithm, at the expense of a small overhead. We 
also mark the cells near hard- wall boundaries as ED cells. Cell bitmasks should be refreshed 
(rebuilt) periodically so as to prevent the fraction of ED cells from increasing. As will be 
seen shortly, it is necessary to introduce at least one "sticky" bit 7^^ that is not cleared but 
rather persists (has memory), and is initialized to zero (not set) at the beginning of the 
simulation. 

2. Near Neighbor Lists 

The cell size should be tailored to the DSMC portion of the algorithm and can become 
much smaller than the size of non-DSMC particles. The LLC method becomes inefficient 
when the interaction (search) range becomes significantly larger than the cell size because 
many cells need to be traversed. In this case the LLC method can be augmented with 
the near-neighbor list (NNL) method, and in particular, the bounding sphere complexes 
(BSCs) method, as described in detail for nonspherical hard particles in Ref. [39]. We have 
implemented the necessary changes to the algorithm to allow the use of NNLs and BSCs (in 
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addition to LLCs), and we used NNLs in our simulations of polymer chains in solution. The 
use of BSCs is not necessary for efficient simulations of polymer solutions if the size of the 
polymer bead is comparable to the size of the cells, which is the case for the simulations we 
report. We do not describe the changes to the algorithm in detail; rather, we only briefly 
mention the essential modiflcations. 

For the purposes of DSMC it is important to maintain accurate particle lists Cc for all cells 
c, so that it is known which particles are in the same cell (and thus candidates for stochastic 
collisions) at any point in time. Therefore, transfers of particles between cells need to be 
predicted and processed even though this is not done in the NNL algorithm described in 
Ref. [39j. Near-neighbor lists are only built and maintained for DSMC particles that are in 
event-driven cells (essentially exactly as described in Ref. [39j). For a DSMC particle i that 
is not in an ED cell q we consider the smallest sphere enclosing cell q to be the (bounding) 
neighborhood (see Ref. |44]) of particle i and only update the (position of the) neighborhood 
when the particle moves to another cell. This ensures that neighbor searches using the NNLs 
are still exact without the overhead of predicting and processing NNL update events for the 
majority of the DSMC particles. 

B. The SEDMD Algorithm 

We have developed an algorithm that combines time-driven DSMC with event-driven MD 
by splitting the particles between ED particles and TD particles. Roughly speaking, only 
the particles inside event-driven cells (i.e., cells for which bit ^ed is set) are part of the AED 
algorithm. The rest of the particles are DSMC particles that are not even inserted into the 
event queue. Instead, they are handled using a time-driven (TD) algorithm very similar to 
that used in classical DSMC. 

It is also possible to implement DSMC as a fully asynchronous event-driven (AED) algo- 
rithm and thus avoid the introduction of an external time scale through the time step At. 
The algorithm introduces a novel type of event we term stochastic (DSMC) collisions^ and 
it is discussed in more detail in Appendix [A] Asynchronous processing has a few advantages 
over the traditional (synchronous) time-driven approach, notably, no errors due to time dis- 
cretization [45j and improved efliciency at low collision rates. For high densities (i.e., high 
collision rates) we have found that these advantages are outweighed by the (implementation 
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and run-time) cost of the increased algorithmic complexity. Additionally, time-driven han- 
dling has certain important advantages in addition to its simplicity, notably, the synchrony 
of the DSMC portion of the algorithm allows for parallelization and easy incorporation of 
algorithmic alternatives (e.g., multi-particle or multi-cell collisions, adaptive open boundary 
conditions, etc.). 

The main types of events in the SEDMD algorithm are: 

Update Move particle i to the current simulation time t if U < t. 

Transfer Move particle i from one cell to another when it crosses the boundary between 
two cells (this may also involve a translation by a multiple of the lattice vectors when 
using periodic BCs). 

Hard-core collision Collide a particle i with a boundary such as a hard wall or another 
particle j with which it interacts. 

Tether collision Bouncing of a pair of tethered particles in a polymer chain when the 
tether stretches (processed exactly like usual hard-particle collisions [38l l40]). 

Time step Move all of the time-driven particles by At and process stochastic collisions 
between them. 

The position r^ and time U as well as the impending event prediction of particle i are updated 
whenever an event involving the particle is processed. 

Both the event-driven and the time-driven DSMC algorithms process stochastic binary 
trial collisions. Processing a trial collision consists of randomly and uniformly selecting a 
pair of DSMC particles i and j that are in the same cell. For hard spheres in the low-density 
limit, the probability of collision for a particular pair ij is proportional to the relative velocity 
Vi^^^ and therefore the pair ij is accepted with probability '^^^^V'^^f'^- If ^ P^i^ accepted 
for collision than the velocities of i and j are updated in a random fashion while preserving 
energy and momentum [30j. If a real collision involving an ED particle i occurs then that 
particle is updated to time txs^ its previous event prediction is invalidated (this may involve 
updating a third-party particle fc), and an immediate update event is scheduled for i (and 
possibly k). 
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It is important to note that the division of the DSMC particles between ED and TD 
handhng is dynamic and does not necessarily correspond to the partitioning of the cells into 
ED and TD cells (based on the cell bitfield ^ed)- As non-DSMC particles move, time-driven 
cells may be masked as event-driven. This does not immediately make the DSMC particles 
in such cells event-driven. Rather, time-driven DSMC particles are moved into the event 
queue only when a collision with a non-DSMC particle is scheduled for them, when they 
move into a TD cell following a time step, or when restarting the event handling. Event- 
driven particles are removed from the event queue when they undergo cell transfer events 
into time-driven cells. 



1. Time Step Events 

The hybrid ED/TD algorithm introduces a new kind of event (not associated with any 
particular particle) called a time step event. This event is scheduled to occur at times 
trs = n/Si, where n G >Z is an integer. When such an event is processed, all of the DSMC 
particles not in the event queue are moved ^ to time trs and are then re-sorted into cells 
(recall that the ED particles are already correctly sorted into cells). Particles that change 
from ED to TD cells and vice- versa are removed or inserted into the event queue accordingly. 
Then, in each cell FcAt trial DSMC collisions are performed, where 

is the DSMC collision rate. Here a — ^"t^R^smc three dimensions and a — ^Rdsmc in 
two dimensions is the coUisional cross-section, Vc is the volume of the cell, and v^ax is an 
upper bound for the maximal particle velocity ^. 

In order to ensure correctness of the AED algorithm, a TD particle must not move by 
more than a certain distance Al^ax when it undertakes a time step. Otherwise, it may 
overlap with a non-DSMC particle that could not have anticipated this and scheduled a 
collision accordingly. Specifically, recall that the event-driven cells are marked whenever a 

^ Note that this update may involve moving some particles by less than At since the time of the last update 
for such particles does not have to be a time step event but could be, for example, a cell transfer. 

^ More precisely, 2vmax is an upper bound on the maximal relative velocity between a pair of particles. 
In our implementation we maintain the maximal encountered particle velocity Vmax and update it after 
every collision and also reset it periodically. 
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neighbor search is performed for a non-DSMC particle. Our simulation uses 

where the masking width Wed '^^ the minimal number of cells covered by any neighbor search 
in any direction, is the (minimal) cell length, and Ddsmc is the diameter of the DSMC 
particles. Any DSMC particle whose velocity exceeds v^ax = ^Imax/^t is inserted into the 
event queue at the end of a time step, and similarly, any particles that would have been 
removed from the event queue are left in the queue if their velocity exceeds the maximum 
safe velocity. Typically, only a small (albeit non-zero) fraction of the DSMC particles falls 
into this category and the majority of the particles that are not in ED cells are not in the 
event queue. In fact, we choose the time step to be as large as possible while still keeping 
the number of dangerously fast DSMC particles negligible. This typically also ensures that 
DSMC particles do not jump over cells from one time step to the next (given that typically 
Wb = 1-2). 

C. Adaptive Open Boundary Conditions 

In three dimensions, a very large number of solvent particles is required to fill the sim- 
ulation domain. The majority of these particles are far from the polymer chain and they 
are unlikely to significantly impact or be impacted by the motion of the polymer chain. It 
therefore seems reasonable to approximate the behavior of the solvent particles sufficiently 
far away from the region of interest with that of a quasi-equilibrium ensemble in which the 
positions of the particles are as in equilibrium and the velocities follow a local Maxwellian 
distribution (the mean of which is equal to the macroscopic local velocity). These particles 
do not need to be simulated explicitly, especially for a DSMC liquid which has no spatial 
structure (ideal gas). Rather, we can think of the polymer chain and the surrounding DSMC 
ffuid as being embedded into an infinite reservoir of DSMC particles which enter and leave 
the simulation domain following the appropriate distributions. 

Such open (Grand Canonical) boundary conditions (BC) are often used in multi-scale 
(coupled) simulations. It is not trivial to implement them when coupling the "reservoir" to 
an MD simulation, especially at higher densities. An example of an algorithm that achieves 
such a coupling for soft-particle systems is USHER [8j. It is also non-trivial to account 
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for the velocity distribution of the particles entering the simulation domain [46], as would 
be needed in a purely event-driven algorithm in which particles are inserted at the surface 
boundary of the domain. However, the combination of a partially time-driven algorithm 
and an unstructured (ideal gas) DSMC fluid makes it very easy to implement open BCs by 
inserting DSMC particles in the cells surrounding the simulation domain only at time-step 
events, based on very simple distributions. 

1. Cell Partitioning 

For the purposes of implementing such non-trivial BCs, we classify the cells as being 
interior, boundary, and external cells. Interior cells are those that are in the vicinity of non- 
DSMC particles, speciflcally, cells that are within a window of half- width Wint > Wed cells 
around the centroid of a non-DSMC particle. The interior cells are divided into event-driven 
and time-driven and are handled as described previously. If a boundary or external cell is 
marked as an event-driven cell the simulation is aborted with an error, ensuring that ED 
cells are always interior. Boundary cells surround the interior cells with a layer of cells of 
thickness Wb > I cells, and they represent cells in which particles may be inserted during 
time step events External cells are non-interior cells that are not explicitly simulated, 
rather, they provide a boundary condition around the interior and boundary cells. This 
layer must be at least Wb cells thick, and the cells within a layer of t^;^ cells around the 
simulation domain (interior together with boundary cells) are marked as both external and 
boundary cells. All of the remaining cells are purely external cells and simply ignored by 
the simulation. Our implementation uses bits in the cell bitmasks to mark a cell as being 
event-driven (bit ^ed)^ boundary (bit 7^), or external (bit 7p). Note that a cell may be a 
combination of these, for example, cells near hard walls might be both interior and boundary, 
and some cells may be both external and boundary. 

Figure [1] provides an illustration of this division of the cells for the simulation of a tethered 
polymer in two and three dimensions. Note that we do not require that the domains of 
interior or non-external cells form a rectangular domain: The flnal shapes and even contiguity 
of such domains depends on the positions of the non-DSMC particles ^. Our implementation 

^ Since ED cells are never boundary cells such insertions cannot lead to overlaps with non-DSMC particles. 
^ If this is not appropriate one can always make the simulation regions (unions of disjoint) rectangular 
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Figure 1: The partitioning of the domain into interior (I) [either event-driven (ED) or time-driven 
(TD)], boundary (B), and external (E) cells in two (left) and three (right) dimensions for a polymer 
chain of 25 beads tethered to a hard wall. The cells are shaded in different shades of gray and 
labeled in the two-dimensional illustration {wed = 2, Wint = 5, wb = 2). The DSMC particles are 
also shown. 

traverses each of the non-DSMC particles in turn and masks the cells in a window of half- 
width w cells around the cell containing the non-DSMC particle as: 

Interior < w < Wint representing cells where the non-trivial flow occurs {wint > wed) 

Boundary Wint < w < Wm + 2wb representing cells where particles may be inserted or 
propagated during time step events. 

External w > Wint + WB representing cells that are not explicitly simulated but rather only 
provide appropriate BCs. 

The division of the cells into event-driven, interior, boundary and external cells is rebuilt 
periodically during the simulation. This rebuilding may only happen at the beginning of 
time steps, and requires a synchronization of all of the particles to the current simulation 
time, a complete rebuilding of the cell bitmasks, and finally, a re-initialization of the event 
processing. Importantly, particles that are in purely external cells are removed from the 
simulation and those that are in event-driven cells are re-inserted into the event queue 

domains simply by padding with interior cells. 
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scheduled for an immediate update event. During the process of rebuilding the cell bitmasks 
cells that are masked as purely external cells are also marked with the sticky bit 75. This 
indicates that these cells need to be re-filled with particles later if they enter the simulation 
domain again (due to the motion of the non-DSMC particles). Once the cell bitmasks are 
rebuilt, a time step event is executed as described next. 

2. Reservoir Particles 

At the beginning of a time step event, after possibly rebuilding the cell masks, the time- 
driven DSMC particles are propagated as usual. If there are external cells, (trial) reservoir 
particles are then inserted into the cells that are both external and boundary, and also in 
cells whose sticky bit 7^ is set (i.e., cells that have not yet been filled with particles), after 
which the bit 7^ is reset. The use of the sticky bit to mark such cells ensures that subsequent 
rebuilding of the masks will not erase the fiag (the sticky bit is only reset once the cell is filled 
with particles). The trial particles are thought to be at a time t — At, and are propagated 
by a time step At to the current simulation time. Only those particles that move into a 
non-external cell are accepted and converted into real particles. If the acceptance would 
insert the particle into a non-boundary cell (i.e., the particle moved by at least vub cells), 
the insertion is rejected and a count of the number of rejected particles reported to aid in 
choosing vub sufficiently large so as to ensure that the tails of the velocity distribution are 
not truncated (in our experience Wb = 2 suffices for reasonable choices of At). Following 
the insertion of reservoir particles stochastic collisions are processed in each cell as usual. 

For improved efficiency, it is possible to replace the volume-based particle reservoir with 
a surface reservoir, and insert particles only at the surface of the simulation domain [46j. 
However, we have not implemented such an approach since the boundary handling is not 
critical for the overall efficiency. 

3. Boundary Conditions 

In our current implementation the reservoir particles follow simple local-equilibrium ideal 
gas distributions. The number of particles to insert in a given cell c is chosen from a Poisson 
distribution with the appropriate density, the positions are uniformly distributed inside 
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the cell, and the velocities are drawn from a biased (local) Maxwellian distribution. The 
mean velocity vm and temperature Tm for the local Maxwellian are chosen according to the 
specified boundary conditions (presently only uniform linear gradients are implemented). 
For example, if a uniform shear in the xy plane is to be applied, vm = TI/c^, where yc is 
the y position of the centroid of the cell and 7 is the shear rate. Using such biased local 
insertions allows one to specify a variety of boundary conditions (for example, a free polymer 
chain in unbounded shear fiow) without resorting to hard-wall boundaries or complicating 
Lee-Edwards conditions. 

It should be noted that in principle we should not use a local Maxwellian velocity distri- 
bution for a system that is not in equilibrium. In particular, for small velocity, temperature, 
and density gradients the Chapman-Enskog distribution is the appropriate one to use in or- 
der to avoid artifacts near the open boundaries at length scales comparable to the mean free 
path A ^J. We judge these effects to be insignificant in our simulations since our bound- 
ary conditions are fixed externally and are thus not affected by the possible small artifacts 
induced in the DSMC fiuid fiow, and since A is small. 

In the future, we plan to replace the particle reservoir with a PDE-based (Navier-Stokes) 
simulation coupled to the DSMC/MD one. Such a fiux-preserving coupling has been im- 
plemented in the past for coupled DSMC/Euler hydrodynamic simulations [471 SH]- It is 
however important for the coupling to also correctly couple fiuctuations. This requires the 
use of fluctuating hydrodynamics in the coupled domain. Such solvers and associated cou- 
pling techniques are only now being developed [26l [27] . 

D. Further Technical Details 

In this section we discuss several technical details of the SEDMD algorithm such as hard- 
wall boundary conditions and the choice of DSMC parameters. 

1. Slip and Stick Boundary Conditions 

We have already discussed the open boundary conditions and their use to specify a variety 
of "far-field" fiow patterns. Additionally, there can also be hard- wall boundaries, i.e., fiat 
impenetrable surfaces. These surfaces can have a velocity of their own and here we discuss 
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how particles reflect from such waUs in the frame that moves with the hard waU. Regardless 
of the details of particle reflections, the total change in linear momentum of all the particles 
colliding with a hard wall can be used to estimate the friction (drag) force acting on the wall 
due to the flow. This can give reliable and quick estimates of the viscosity of a DSMC fluid, 
for example. We use the classical no-slip BCs (i.e., zero normal and parallel velocity) for 
smooth hard- wall surfaces. Molecular simulations have found some slip; however, at length- 
scales signiflcantly larger than the mean free path and/or the typical surface roughness one 
may assume no-slip boundaries if the hard-wall boundary position is corrected by a slip 
length Lsiip [3j. 

Our simulations of tethered polymers use thermal walls (kept at kT = 1) [30j to implement 
no-slip hard walls at the boundaries of the simulation cell. Following the collision of a particle 
with such a wall, the particle velocity is completely randomized and drawn from a half 
Maxwell-Boltzmann distribution (other biased distributions may be used as appropriate). 
This automatically ensures a zero mean velocity at the wall boundary and also acts as a 
thermostat keeping the temperature constant even in the presence of shear heating. No-slip 
boundaries can also be implemented using (athermal) rough walls which reflect incoming 
particles with velocity that is the exact opposite of the incoming velocity [49j . Similarly, slip 
boundary conditions (zero normal velocity) can be trivially implemented by using specular 
walls that only reverse the normal component of the velocity (relative to the wall). A 
mixture of the two can be used to implement partially rough walls, for example, a roughness 
parameter < < 1 can be used as the probability of randomly selecting a rough versus 
a specular collision. 

Similar considerations apply to the boundary conditions at the interface of a hard particle 
such as a polymer bead. Most particle-based methods developed for the simulation of particle 
suspensions consider the solvent particles as point particles for simplicity, and only MD 
or certain boundary discretization schemes [50] resolve the actual solvent-solute interface. 
Specular BCs are typical of MD simulations and assume perfectly conservative collisions 
(i.e., both linear momentum and energy are conserved). However, if the polymer beads are 
themselves composed of many atoms, they will act as a partially thermal (and rough) wall 
and energy will not be conserved exactly. 

In the simulations reported here we have used rough walls for collisions between DSMC 
and non-DSMC particles. This emulates a non-stick boundary condition at the surface of 
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the polymer beads. Using specular (slip) conditions lowers the friction coefficient ^, but does 
not appear to qualitatively affect the behavior of tethered polymers. 

2. Constant Pressure Flows 

We note briefly on our implementation of constant pressure boundary conditions, as used 
to simulate flow through open pipes. A constant pressure is typically emulated in particle 
simulations via a constant acceleration a for the DSMC particles [51j together with pe- 
riodic BCs along the flow (acceleration) direction. In time-driven algorithms, one simply 
increments the velocity of every particle by a At and the position by aAt^/2 at each time- 
step (before processing DSMC collisions). In SEDMD this is not easily implemented, since 
the trajectory of the DSMC particles becomes parabolic instead of linear and exact colli- 
sion prediction between the DSMC and the non-DSMC particles is complicated. We have 
opted to implement constant pressure BCs by using a periodic delta-function forcing on the 
DSMC particles. Speciflcally, the velocities of all DSMC particles are incremented at the 
beginning of each time step by aAt, and then stochastic collisions are processed. 

3. Choice of DSMC Collision Frequency 

The viscosity of the DSMC fluid is determined by the choice of collision frequency Tc 
and cell size Lc. Classical DSMC wisdom [30j is that cell size should be smaller than the 
mean free path, Lc <C A, but large enough to contain on the order of Nc ^ 20 particles 
(in three dimensions). It is obvious that both of these conditions cannot be satisfled for 
denser liquids, where A is only a fraction of the particle size. It is now well-known that it 
is not necessary to have many particles per cell, so long as in Eq. ([T]) we use Nc{Nc — 1) 
instead of the traditional (but wrong) N^. Coupled with the Poisson distribution of Nc this 
gives a constant average total collision rate. However, using very small cells leads to very 
large variability of collision rates from cell to cell and thus spatial localization of momentum 
transfer during each time step. Namely, with very small cells one rarely has two particles and 

^ The Stokes friction force has a coefficient of 47r for slip BCs instead of the well-known Gtt for no-slip BCs. 
Recall that the event prediction for any ED particle i whose velocity is changed must be updated, typically 
by scheduling an immediate update event. 
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thus most of the coUisions wiU occur in the few ceUs that happen to be densely populated. 

We have aimed at trying to mimic what would happen in an MD simulation in the DSMC 
one. In an MD simulation particles collide if their distance is equal to the particle diameter 
D. Therefore, we have aimed at keeping the cell size at a couple of diameters, Lc ^ 2D. 
At typical hard-sphere liquid densities this leads to A^c ^ 5 — 10, which seems appropriate 
in that it allows enough collision partners for most of the particles but still localized the 
momentum transfer sufficiently. For very small mean free paths DSMC does not distinguish 
velocity gradients at length scales smaller than the cell size and, in a long-time average sense, 
localizes the velocity gradients at cell interfaces We will assume that, for problems of 
interest to us, the structure of the fluid and flow at length-scales comparable to D (and thus 
Lc) is unimportant, and verify this by explicit comparisons to MD. 

When the cell size is chosen such that Nc ^ 5 — 10 and the time step is reasonable. 
At ^ (0.1 - 0.2)Ljv, Eq. Q gives collision frequencies that are sufliciently high so that 
almost all particles suflFer at least one collision every time step, and typically more than 
one collision. The effect of such repeated collisions is to completely thermalize the flow to 
a local equilibrium (i.e., local Maxwellian), and we have observed that further increasing 
the collision frequency does not change the effective viscosity (we do not have a theoretical 
understanding of this behavior t52j). For greater efl&ciency, we have chosen to use the lowest 
collision rate (for a given timestep) that still achieves a viscosity that is as high as using a 
very high collision rate. We flnd that this is typically achieved when each particle suffers 
about half a collision or one collision each timestep [53] . Appendix [B] describes some multi- 
particle collision variants that may be more appropriate under different conditions. 

4. DSMC without Hydrodynamics 

The solvent exerts three primary effects on polymers in flow: (1) stochastic forces due to 
fluctuations in the fluid (leading to Brownian-like motion), (2) (local) frictional resistance 
to bead motion (usually assumed to follow Stokes law), and (3) hydrodynamic interactions 
between the beads due to perturbations of the flow fleld by the motion of the beads. Brow- 
nian dynamics, the most common method for simulating the behavior of polymers in flow, 
coordinates the flrst two effects via the fluctuation-dissipation theorem, and potentially adds 
the third one via approximations based on the Oseen tensor (neglecting the possibility of 
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large changes to the flow fleld due to the moving beads). By turning oS local momentum 
conservation one can eliminate all hydrodynamic interactions, and thus test the importance 
of the coupling between polymer motion and flow. 

Yeoman's et ai [20l [34] have implemented a no-hydrodynamics variant of the MPCD 
algorithm by randomly exchanging the velocities between all particles at each time step (thus 
preserving momentum and energy globally, but not locally). In the presence of a background 
flow, such as shear, only the components of the velocities relative to the background flow 
are exchanged. We have implemented a no-hydrodynamics variant of DSMC by neglecting 
momentum conservation in the usual stochastic binary collisions Speciflcally, if particles 
A and B collide, the post-collisional velocity of A is set to be the same magnitude as that of B 
but with a random orientation, and vice versa (this conserves energy but not momentum). 
If the boundary conditions specify a background flow such as a uniform shear the flow 
velocity is evaluated at the center of the DSMC cell and the collisions are performed in the 
frame moving with that velocity. This forces the average velocity proflle to be as specifled 
by the boundary conditions, but does not allow for perturbations to that proflle due to 
hydrodynamic effects. 

IV. PERFORMANCE IMPROVEMENT 

It is, of course, expected that the DSMC algorithm will give a performance improvement 
over MD. However, to make an impact on real- world problems this performance gain must be 
an order of magnitude or more improvement. Indeed, we flnd that SEDMD with adaptive 
boundary conditions can be up to two hundred times faster than EDMD under certain 
conditions. Note also that it is well-known that EDMD is already signiflcantly faster than 
TDMD (depending on the density), although such a comparison is somewhat unfair since the 
hard-core interaction potentials are very simple by design. It is important to note that this 
algorithm is serial, and we do not consider or use any parallelization. Because of the inherent 
simplicity and thus efliciency of the algorithm, however, it is possible to study time scales 
and system sizes as large or larger than parallel simulations described in the literature so 
far. The combined time-driven DSMC with event-driven MD algorithm can be parallelized 

In this implementation switching hydrodynamics off becomes an alternative branch localized in the binary 
collision routine and the algorithm is otherwise unchanged. 
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using traditional techniques from TDMD if proper domain partitioning is constructed, so 
that each event-driven region is processed by a single processor 

As model problem we study a tethered polymer in three dimensions. The solvent density 
was chosen to be typical of a moderately dense hard-sphere liquid. The performance and 
optimal choice of parameters depends heavily on the size of the beads relative to the size 
of the solvent particles for both MD and the hybrid algorithm. Realistically, beads (meant 
to represent a Kuhn segment) should be larger than the solvent molecules^^. This of course 
dramatically increases the computational requirements due to the increase in the number 
of solvent particles (and also makes neighbor searching more costly). For this reason, most 
MD simulations reported in the literature use solvent particles that are equivalent, except 
for the chain connectivity, to the solute particles. 

Our first test problem is for a chain of 25 large beads, each about 10 times larger (in 
volume and in mass) than the solvent particles, in a box of size 2 x 1.25 x 1.25 polymer 
lengths, for a total of about = 2.3 x 10^ particles. For the DSMC simulations, we did not 
use BSCs (bounding sphere complexes [39]), and therefore the neighbor search had to include 
next-nearest neighbor cells as well (i.e., Wed = 2). For the corresponding MD simulations, 
BSCs were used. Under these conditions, DSMC outperformed MD by a factor of 35. If 
adaptive open BCs were used with Wm = 5, giving about N = 3.2 x 10^ particles (the exact 
number changes with polymer conformation), the speedup was 180. While this may seem an 
unfair comparison, it is important to point out that we do not even know how to implement 
an adaptive simulation domain in pure EDMD. 

The second test problem was for a chain of 30 beads which were identical to the solvent 
particles, except for the added chain tethers. The number of particles in the simulation cell 
was thus much smaller, N = 4.8 x 10^, and Wed = 1- Adaptive BCs with Wm = 5 reduce the 
simulation domain to A = 2.2 x 10^ particles. For these parameters DSMC with adaptive 
BCs was about 30 times faster than full MD. Table |T] summarizes the large performance 
gains of SEDMD relative to traditional EDMD. 

One of the fundamental problems with multi-scale modeling is that typically the ma- 
jority of the simulation time is spent in the finest model since it is difficult to match the 

Achieving good load balancing will be easiest for systems containing multiple polymer chains. 

For example, in Ref. [19] an appropriate bead size for polyethylene is estimated at 1.5nm, and for DNA 

(a much stiffer molecule with large persistence length) at 40nm. 
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Standard BCs 


Adaptive BCs 


Large beads 


35 


180 


Small beads 


20 


30 



Table I: Performance gains of SEDMD relative to EDMD for a typical tethered polymer simulation. 

time scales of the coupled components [24j. For example, MD simulations are so expensive 
that coupling them to almost any meso- or macro-scopic solver leads to simulation times 
limited by that of MD simulations (albeit of a much smaller system). By virtue of the fast 
microscopic algorithm (EDMD instead of TDMD) and the efficient coupling, our method 
spends comparable amounts of computation on the solute (and immediately surrounding 
solvent) and the solvent particles. For the DSMC run with adaptive open BCs and large 
beads, about 50% of the time was spent in manipulation of near neighbor lists. Most of 
the remaining time was spent inside the routine that takes a DSMC timestep, and actual 
processing of DSMC collisions (both trial and real) occupied about 20% of the computation 
time. For small beads, the majority of the time, 80%, was spent in the DSMC time-step 
routine, and processing of DSMC binary collisions occupied about 35% of the computation 
time. 

V. TETHERED POLYMER IN SHEAR FLOW 

In this section results are presented for a tethered polymer chain in uniform shear in three 
dimensions. The linear chain is in a good solvent and is attached at one end to a hard wall, 
as represented by the plane y = 0. A linear velocity profile v = ^yx along the x axis is 
imposed sufficiently far from the chain. This problem was first studied experimentally by 
Doyle et al [7J and since then numerous computational studies have investigated various 
aspects of the problem [HI |9l [101 [HI [IS] . We will focus on the dynamics of the chain at low 
to medium ffow rates (Weissenberg numbers) because we wanted to verify that our polymer 
and solvent model can correctly reproduce non-trivial dynamics. 
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A. Background 

The properties of a linear polymer in shear flow can be related to the dimensionless 
Weissenberg number Wi = 7ro, where tq = r(7 = 0) is the relaxation time of the polymer 
chain when there is no shear. When Wi < 1 the flow barely affects the polymer, contrary 
to when Wi > 1. Different models have given similar properties for the same Weissenberg 
number. 

The original experimental study of tethered polymers [7] observed what was termed "cyclic 
dynamics" of the chains. Speciflcally, the following cycle was proposed. When the polymer 
moves too far from the wall, presumably by an unusual fluctuation, it experiences a stronger 
flow and is stretched. A torque develops that then pushes the chain closer to the wall, 
where it can contract again due to the weaker flow near the wall. The cycle then repeats. 
Experiments [7J did not identify clear periodicity of this motion. Subsequent computational 
studies have looked for such a characteristic period for this cycling motion. 

The MD study in Ref. [9j examined the cross-correlation function Cx(f){t)^ where X 
measures the extension of the polymer along the flow, and measures the angle of the chain 
with respect to the hard wall. No exact deflnitions of A or were given even though there 
are several possibilities. One can use the difference between the maximal and the minimal 
bead positions as a measure of the extension along a given axes. Optionally, one can simply 
use the maximal position, or one can use the position of the last bead. Similarly, the angle 
of the polymer can be based on a linear flt to the shape of the chain, on the position of the 
center of mass, the asymmetry of the gyration tensor [12j, or the position of the last bead. 
We have examined various choices and have found little qualitative difference between the 
different choices. We have found the position of the end bead r^y^ = (x, ^, z) to be the best 
option and will also measure the angle = tan~-^(y/x). 

The authors of Ref. [9j found that Cx(f){t) develops a peak at positive time t* for sufficiently 
large Wi numbers (Wi > 10). This was interpreted as supporting the existence of a critical 
Weissenberg number Wi where the flow effect on the polymer dynamics changes qualitatively. 
It was also found that t* decreases with increasing Wi and the height of the peak increases. It 
is important to note that t* was found to be comparable to the relaxation time of the polymer 
tq. Additionally, the internal relaxation time r was found to decrease with increasing Wi, 
in agreement with theoretical predictions. 
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A subsequent study which used a hybrid MD/CFD model, and also a (free-draining) 
Brownian dynamics model, claimed to observe periodic oscillations in the cross-correlation 
function between the extensions along the flow and along the shear direction (i.e., perpen- 
dicular to the wall), Cxy{t) [lOl [13]. However, the period of oscillation was found to be an 
order of magnitude larger than the internal relaxation time, as revealed by a small peak in 
the power spectral density PSDxy{f) of Cxy{t). A similar claim was made in Ref. [12j based 
on PSDcfycfy of the polymer angle autocorrelation function Ccf)cf){t) for both a free polymer 
in unbounded shear flow and a tethered polymer in shear flow. No results for the short- 
time cross-correlation functions were reported in either of these studies making it diflicult 
to reconcile the results obtained from PSDs with those in Ref. [9j . 

Most experimental and computational studies of the dynamics of polymers in shear flow 
have been for free chains in unbounded flow [4j. In that problem, for Wi > 1, it is possible 
to identify a well-deflned "tumbling" event as the polymer rotates. The frequency of such 
tumbling times can be measured by visual inspection and have been compared to the com- 
puted location of the peak in the PSDs ^2l[55- The good match has thus been taken as an 
indicator that PSDs peaks can be used to determine characteristic tumbling times and the 
same methodology has been applied to a tethered polymer as well. However, for the case of 
a tethered chain it is not easy to identify a periodic event such as a speciflc rare fluctuation. 
Therefore, it is not surprising that we do not conflrm the existence of a characteristic time 
that is an order of magnitude larger than the internal relaxation time. One must here dis- 
tinguish between "cyclic" (repetitive) events and periodic events. A Poisson time process of 
rate P has a well-deflned time scale P"^, however, the occurrence of such events is not peri- 
odic; the delay between successive events is exponentially-distributed. In Ref. [54j such an 
exponential distribution is proposed even for the delay between successive tumbling events 
for a free chain in unbounded flow. The PSD of such a process is expected to be that of 
white noise (i.e., flat) for frequencies small compared to P, and typically a power-law decay 
for larger frequencies (gray noise). The occurrence and shape of any local maxima (peaks) 
or frequencies comparable to P depends on the exact nature of the correlations at that time 
scale. 

The PSD is equivalent to the Fourier spectrum power of the angle trace (j){t) based on the convolution 
theorem. 
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B. Model Parameters 



As explained in [TV] we have made several runs for different polymer lengths and also 
bead sizes. One set of runs used either A^^ = 25 or 50 large beads each about 10 times larger 



than a solvent particle, using DSMC with or without hydrodynamics (see Section HID 4) 
for the solvent. Another set of runs used either A5 = 30 or 60 small beads each identical to a 
solvent particle, using DSMC or pure MD for the solvent. The beads were rough in the sense 



that no-slip conditions were applied for the solvent-solute interface (see Section HID 1). All 



of the runs used open boundary conditions (see Section IIIC), and the typical half- width of 
the interior region was Wm = 5 or Wint = 7 cells around the polymer chain. The difference 
in the results (such as relaxation times) between these runs and runs using = 10 or 
runs using periodic BCs were negligible for the chain sizes we studied The solvent was a 
hard-sphere MD or DSMC fluid with volume fraction (j) ^ 0.25 — 0.30, which corresponds to 
a moderately dense liquid (the melting point is 0^ ^ 0.49). The A5 = 30 runs were run for 
T ^ 6000ro with Wint = 7, and such a run takes about 6 days on a single 2.4GHz Dual-Core 
AMD Opteron processor. Even for such long runs the statistical errors due to the strong 
fluctuations in the polymer conformations are large, especially for correlation functions at 
long time lags t > r. 



C. Relaxation Times 



The relaxation time of the polymer r is well-defined only for linear models. It is often 
measured by fitting an exponential to the autocorrelation function of the end-to-end vector 
^end{t) = r^v^ — ri, where denotes the position of the i-th bead [6j. We will separately 
consider the different components of the end-to-end vector r^nd = {^^Vi^) and fit an expo- 
nential to the Cxx^ Cyy and Czz auto-correlations functions to obtain the relaxation times 

It is expected that using a small Wint would truncate the (long-ranged) hydrodynamic interactions and 
thus increase the relaxation time. We observe such effects for the N^j = 50 chains, however, the effect is 
too small compared to the statistical errors to be accurately quantified. 

The initial relaxation of the various auto-correlation functions C{t) is faster than exponential, and the 
statistical error at longer times is large even for long runs. We therefore fit the exponentials to the portion 
of the curves at small times, when 0.2 < C{t) < 0.8. The fits are not perfect and there are large statistical 
errors depending on the length of the run and the number of samples used to average C(t), and the 
relaxation times (and thus Weissenberg numbers) we quote should be taken as approximate. 
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Figure 2: Dependence of the relaxation times of the different components of the end-to-end dis- 
placement vector on the Weissenberg number. The relaxation times have been renormalized to 
equal unity for Wi = for direct comparison. For each Wi, Tx is shown with circles, Ty with 
squares, and with diamonds. Different textures of the symbols are used for the different models, 
as indicated in the legend. The inset shows Ty/rx and Tz/tx for the different runs. 

Tx^ Ty and Tz as a function of Wi. We find that is always the largest, especially for large 
Wi (for Wi = 0, = Tx by symmetry), and Ty is always smaller by at least a factor of two ^ 



even for Wi = 0, as illustrated in the inset in Fig. [2] We take tq = Tx{Wi = 0) = r^(Wi = 0) 
as the definition of the polymer relaxation time. 

The relaxation times we observe for Wi = are consistent with what is predicted from 
theoretical considerations, r ^ 0.9r]b^N^-^ /kT^ where r] is the viscosity and b is the effective 
bead radius. Using the viscosity (based on Enskog theory) of the MD liquid and the tether 
length as 6, we calculated t ^ 19 for the case of A^^ = 25 with large beads, to be compared 



This is because of the constraint that the polymer chain must be above the plane ^ = at all times, which 
reduces the available configuration space. 

Direct measurements of the viscosity of the DSMC liquid show that it has viscosity rather close to that 
of the corresponding MD liquid for the specific parameters we use. 
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to the numerical results from DSMC r = 25±5. The MD runs for the case of large beads are 
not long enough to determine the relaxation time accurately. We expect that the difference 
between MD and DSMC will become more pronounced for smaller beads, and indeed, for 
= 30 we obtain tmd ^ ^^dsmc- Turning hydrodynamics off in DSMC extends the 
relaxation times (and also the collapse times for an initially stretched polymer) by a factor 
of 3 — 5, as already observed using MPCD [20j and as predicted by Zimm theory Figure 
[2] illustrates the dependence of Tx{Wi)/Tx{Wi = 0) on Wi, and similarly for the y and z 
directions. Quantitatively similar (but not identical) results are observed independently of 
the details of the polymer model and even the existence of hydrodynamic relaxations. 

D. Cyclic Dynamics 

We now turn our attention to cross-correlations between polymer extensions in the x and 
y directions. We have found that the cross-correlations lags are most visible in the x and 
y positions of the last bead, Cxy{t). Our results for Cxy{t) are shown in Fig. |3| along with 
Cx(j){t) as an inset. The results for Cx(f){t) compare well with those in Ref. [9j, although we 
see the secondary peak developing at somewhat lower Wi. We do not see any evidence for 
the existence of a critical Wi: There are peaks at both positive and negative time in Cxy{t) 
for all Wi. Some cross-correlations, such as Cx(f){t)^ have a large positive or negative cusp at 
the origin at Wi = and it is this cusp that masks the peaks at non-zero lags for small Wi. 

In Fig. [4] we compare Cxy{t) at Wi ^ 2 for several different models and see a good 
match, even for the DSMC runs ignoring hydrodynamics (momentum conservation). This 
indicates that the dynamics of the chains is primarily driven by the competition between 
the internal stochastic motion (entropy) and the external forcing due to the shear, and not 
hydrodynamic interactions between the beads or the effect of the motion of the chain on the 
flow. 

We do not discuss the origin and locations of the peaks in the cross-correlation functions 
in detail in this work. These peaks are indicative of the existence of a correlated motion in 

It is difficult to directly compare DSMC with and without hydrodynamics since switching hydrodynamics 
off, in our model, affects the friction force between the beads and the solvent. This is unhke the models 
were the friction force is an added phenomenological term that has an adjustable coefficient. 
The Weissenberg numbers were calculated after the runs were completed and therefore the different runs 
are not at the exact same Wi number. 
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Figure 3: Cross correlation function Cxy{t) for chains of N^, = 30 small beads in a DSMC solvent 
at different shear rates. The inset shows the corresponding Cx(f)(t) for comparison with the soft- 
particle MD results in Ref. [9j. Peaks are visible in Cxy{t) at all Wi > 0, but are obscured in 
Cx(j){t) due to the large negative dip at the origin for Wi = 0. There are large statistical errors at 
small Wi making it difficult to identify the peaks. 

the xy plane, but do not uniquely identify that motion. An important question to address 
is the existence of a time scale other than the internal relaxation time r(Wi). In Fig. [5] we 
show a renormalized cross-correlation function 



a 



xy 



Wi 



t 



_r(Wi)_ 

in an unsuccessful attempt to collapse the data for different Wi. While the match is not 
perfect the picture does not point to the existence of a time scale shorter than r(Wi). We 
also do not see any convincing evidence for coherent and reproducible correlations on time 
scales significantly larger than r, even in various power spectral densities. Our results do not 
rule out the possibility of a repetitive motion of the chain with widely varying cyclic times 
(e.g., exponential tail) but we have not observed any direct evidence for such cycling either. 
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Figure 4: Comparison of Cxy{t) for Weissenberg number of about 2 for several different models, 
after the time axes has been normalized. 

We will report more detailed results on the dynamics of tethered polymer chains along with 
comparisons with Brownian dynamics and Lattice-Boltzmann in future work. 

VI. CONCLUSIONS 

We presented a stochastic event-driven molecular dynamics (SEDMD) algorithm that 
combines hard-sphere event-driven molecular dynamics (EDMD) with direct simulation 
Monte Carlo (DSMC), aimed at simulating flow in suspensions at the microscale. The 
overall algorithm is still event-driven, however, the DSMC portion of the algorithm can be 
made time-driven for increased efliciency. The fundamental idea is to replace the determin- 
istic (MD-like) interactions between particles of certain species with a stochastic (MC-like) 
collision process, thus preserving the phase space dynamics and conservation laws but ig- 
noring the liquid structure. The SEDMD methodology correctly reproduces hydrodynamic 
behavior at the macroscale but also correctly represents fluctuations at the microscale. A 
similar algorithm has been proposed using time-driven (soft-particle) MD and a multiparticle 
collision variant of DSMC [23] . 

As an application of such a methodology we have considered the simulation of polymer 
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Figure 5: Cross correlation function Cxy for = 30 DSMC runs as in Fig. [s] but with time 
renormalized by r{Wi) and the correlation magnitude scaled by Wi. 

chains in a flowing solution, and in particular, a polymer tethered to a hard wall and subject 
to shear flow. We have implemented open boundary conditions that adaptively adjust the 
simulation domain to only focus on the region close to the polymer chain(s). The algorithm 
is found to be eflicient even though it is not parallelized, and it is found to reproduce results 
obtained via molecular dynamics and other algorithms in the literature, after adjusting for 
the correction to transport coeflicients and compressibility of the DSMC fluid relative to the 
MD fluid. 

We studied the dynamics of a tethered polymer subject to pure shear and found consistent 
results between MD and DSMC and also previous TDMD studies. We flnd that neither the 
size of the polymer beads relative to the solvent particles, nor the correct representation of 
the hydrodynamic interactions in the fluid, qualitatively alter the results. This suggests that 
fluctuations dominate the dynamic behavior of tethered polymers, consistent with previous 
studies. Our results do not flnd periodic motion of the polymer and show that the cross- 
correlation between the polymer extensions along the flow and shear directions shows a 
double-peak structure with characteristic time that is comparable to the relaxation time of 
the polymer. This is in contrast to other works that claim the existence of a new timescale 
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associated with the cychc motion of the polymer. We wiU investigate these issues further 
and compare with Brownian dynamics and Lattice-Boltzmann simulations in future work. 

We expect that this and related algorithms will find many applications in micro- and 
nano-fiuidics. In particular, the use of DSMC instead of expensive MD is suitable for prob- 
lems where the detailed structure and chemical specificity of the solvent do not matter, and 
more general hydrodynamic forces and internal fiuct nations dominate. Using a continuum 
approach such as Navier-Stokes (NS) equations for the solvent is questionable at very small 
length scales. Furthermore, the handling of singularities and fiuctuations is not natural in 
such PDE methods and various approximations need to be evaluated using particle-based 
methods. Since the meshes required by continuum solvers for microfiows are very fine, it is 
expected that the efficiency of particle methods will be comparable to PDE solvers. Nev- 
ertheless, algorithms based on ffuctuating hydrodynamics descriptions will be more efficient 
when ffuctuations matter. Comparisons and coupling of DSMC to ffuctuating NS solvers is 
the subject of current investigations [27j. 
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Appendix A: AED VARIANTS OF DSMC 

In this Appendix we discuss a fully asynchronous event-driven (AED) implementation of 
DSMC. The advantage of asynchronous algorithms is that they do not introduce any artifi- 
cial time scales (such as a time step) into the problem [55j . We have validated that the AED 
algorithm produces the same results as the time-driven one by comparing against published 
DSMC results for plane Poiseuille fiow of a rare gas [51j- We have also implemented tradi- 
tional time-driven (TD) DSMC and find identical results when the time step is sufficiently 
small. We find that the event-driven algorithm is almost an order of magnitude slower than 
the time driven one at higher densities, and only becomes competitive at very low densities 
(which is the traditional domain of interest for DSMC). The overhead of the AED algorithm 
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comes from the need to re-predict the next event and update the event queue whenever 
a particle suffers a DSMC coUision. This cost is in addition to the equivalent cost in the 
time-driven algorithm, namely, moving the particles forward in time and colliding them. 
The AED algorithm introduces a new type of event, a stochastic (trial) collision between 



two DSMC particles that are in the same cell (see Section IIIB). These trial collisions occur 
in a given cell c as a Poisson process (i.e., exponentially distributed waiting time) with a 
rate given by Eq. ([T]). There are several approaches to scheduling and processing DSMC 
collisions directly borrowed from algorithms for performing Kinetic Monte Carlo simulations 
(which are synchronous event-driven algorithms ^6]). The simplest, and in our experience, 
most efficient, approach to AED DSMC is to use cell rejection to select a host cell for 
the stochastic collisions. The rate of DSMC collisions is chosen according to the cell with 
maximal occupancy N^^^^ T = Nceiis^T^^ • The randomly chosen cell c of occupancy Nc 
is accepted with probability Nc{Nc — 1)/ ]^N^^^{N^^^ — 1)] and a random pair of particles 
i and j are chosen from Cc- Since the DSMC fluid is perfectly compressible, the maximal 
cell occupancy can be quite high for very large systems, and this leads to decreasing cell 
acceptance probability as the size of the system increases. 

One can avoid cell rejections altogether. The first option is to associate stochastic colli- 
sions with cells and schedule one such event per cell. The event time is easily predicted at 
any point in time t to occur at time t — F"^ Inr , where r is a uniform random deviate in 
(0, 1). These event times are put in the event queue (which may be separate from the particle 
one and then the two queues may be merged only at the top). The event times (and thus 
the cell queue) need to updated whenever a cell occupancy Nc changes, that is, whenever 
a cell transfer is processed. This makes this algorithm inefficient. Another alternative is to 
recognize that the sum of a set of independent Poisson processes is a Poisson process with 
a rate that is the sum of the individual rates, P = "^^c^c- That is, DSMC collisions occur 
in the system as a Poisson process with rate P. When processing such an event one has to 
first choose the cell with probability Pc/P, which requires some additional data structures 
to implement efficiently ^6]. For example, the cells could be grouped in lists based on their 
occupancy and then an occupancy chosen first (with the appropriate weight), followed by 
selection of a cell with that particular occupancy. 

Finally, it is also possible to use a mixture of the asynchronous and time-driven variants 
of DSMC. The asynchronous algorithm (e.g., based on cell rejection) can be used for DSMC 
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particles in event-driven cells, and the time-driven one elsewhere. This may be useful in 
situations where the time-scale of the event-driven component (i.e., the solute and nearby 
solvent particles) is significantly smaller than the time step At and thus time stepping would 
lead to discretization artifacts. 



In the AED variant of DSMC constant pressure BCs (see Section HID 2) can be im- 
plemented by adding a new type of acceleration event. When such an event is processed, 
all of the particles are brought to the same point in time (synchronized), the velocities of 
each DSMC particle i is incremented by aAt^ (here Ati is the elapsed time since the last 
acceleration event), and the event handling is restarted. The acceleration events occur as a 
Poisson process with a suitably chosen rate, for example, ensuring that the average or max- 
imal change in velocity is a fraction of the average particle velocity. Note that the choice of 
this acceleration rate introduces an artificial time constant in the algorithm similar to the 
time step At in time-driven DSMC. 

Appendix B: MULTI-PARTICLE COLLISIONS IN DSMC 

Under dense liquid conditions, DSMC binary collisions are so numerous (see Section 



III D 3 ) that the velocities of the particles are effectively thermalized to the local Maxwell 
distribution. We have implemented a variant DSMC algorithm in which at every time step 
the velocities of all of the particles are redrawn from a local Maxwellian, preserving the 
total linear momentum and energy in each cell [57j. We found that this variant of DSMC is 
less efficient than and behaves similarly to the usual binary-collision DSMC. Reference [58] 
describes a more general algorithm (TRMC) that combines binary collisions for a subset 
of the particles with drawing from a local Maxwellian for the remainder of the particles, 
and under dense liquid conditions this typically degenerates to complete randomization of 
all of the velocities at every time step. Until a theoretical framework is established for the 
behavior of DSMC-like algorithms at high densities the classical DSMC algorithm seems to 
be the best alternative in terms of simplicity, efficiency, and theoretical foundation. The 
effect of collision rules and cell size on multiparticle collision dynamics has been studied and 
it was found that increased collisional viscosity is desirable for achieving realistic convection 
to diffusion ratios [21 [22] . 

We mention that, strictly speaking, we should use as K in Eq. ([T]) not the volume of 
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the cell, but the unoccupied cell volume (that is, the portion of the cell not covered by non- 
DSMC particles) It is however difficult to dynamically maintain an accurate estimate 
of the cell coverage, and the complication does not appear to be worth the implementation 
complexity. In particular, an approximation is already made in neglecting the structure of 
the fluid near a polymer bead (i.e., the solvation layer), and furthermore, the majority of the 
cells that are partially covered by a polymer bead will be entirely or almost entirely covered 
so that they would at most contain a single DSMC particle, in which case the probability 
of a DSMC collision would be very low anyway. Finally, as explained in Section |IIID 3| the 
exact collision frequency does not really matter. In the context of multiparticle collision 
dynamics, Ref. [49j proposes the use of virtual particles fllling the partially-flUed cells as a 
way to achieve more accurate stick boundary conditions. 
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